home *** CD-ROM | disk | FTP | other *** search
/ Programmer Power Tools / Programmer Power Tools.iso / progjrn / pj_6_6.arc / CCODE.ARC / FLTBENCH.C < prev   
Encoding:
C/C++ Source or Header  |  1988-09-01  |  37.7 KB  |  1,465 lines
  1. #ifdef __HIGHC__
  2. #pragma wtrsvd1(16);    /* New Intel setting for Weitek registers. */
  3. #endif
  4. /* ------------------------------------------------------------------
  5.  
  6.     Floating-point benchmarks.
  7.     
  8. -------------------------------------------------------------------*/
  9.  
  10.  
  11. #include <stdio.h>
  12. #include <math.h>
  13. #include "bench.h"
  14.  
  15. #define TRUE     1
  16. #define FALSE    0
  17.  
  18. /*
  19.  
  20. main - Main driver module for benchmarks
  21.  
  22. */
  23.  
  24. null_init(){}
  25.  
  26. main(L_argc,L_argv)
  27. int L_argc; char ** L_argv;
  28.  
  29. {
  30.  
  31.     extern Float();
  32.     extern flt();
  33.     extern savage();
  34.     extern autodesk(), autodesk_init(), INTRIG_autodesk();
  35.     extern whetd(),whets();
  36.     long Wtime;
  37.  
  38.         argc = L_argc; argv = L_argv;
  39.     printf("\n");
  40.     printf("Floating-point benchmarks from BYTE Magazine (July 1987, Page 101),\n");
  41.     printf("and other standard benchmarks -- version 1.1, corrected Whetstones.\n");
  42.     printf("\n");
  43.     common_header();    /* Dobench knows the format of this. */
  44.     DoBench("F", "Float", 10000L, Float, null_init, TRUE);
  45.     DoBench("f", "Flt", 256000L, flt, null_init, TRUE);
  46.     DoBench("S", "Savage", 25000L, savage, null_init, TRUE);
  47.         DoBench("A", "Autodesk",1000L, autodesk, autodesk_init, TRUE);
  48.         DoBench("AI", "Autodesk (INTRIG)",1000L, INTRIG_autodesk, autodesk_init, TRUE);
  49.         Wtime = DoBench("WS", "Whetstone (single)", 100L, whets, null_init, FALSE);
  50.            if (Wtime) printf("  (%d KWhets/sec).\n",(int)(10000 / (Wtime/100.0))); 
  51.         Wtime = DoBench("WD", "Whetstone (double)", 100L, whetd, null_init, FALSE);
  52.            if (Wtime) printf("  (%d KWhets/sec).\n",(int)(10000 / (Wtime/100.0)));
  53.            
  54.     common_trailer();
  55. }
  56.  
  57. /*------------------------------------------------------------------*/
  58.  
  59. /* simple benchmark for testing floating point speed of c libraries
  60.    does repeated multiplications and divisions in a loop that is
  61.    large enough to make the looping time insignificant */
  62.  
  63. #define CONST1  3.141597E0
  64. #define CONST2 1.7839032E4
  65. #define COUNT 10000
  66.  
  67. Float()
  68.         {
  69.         double a, b, c;
  70.         int i;
  71.  
  72.         a = CONST1;
  73.         b = CONST2;
  74.         for (i = 0; i < COUNT; ++i)
  75.                 {
  76.                 c = a * b;
  77.                 c = c / a;
  78.                 c = a * b;
  79.                 c = c / a;
  80.                 c = a * b;
  81.                 c = c / a;
  82.                 c = a * b;
  83.                 c = c / a;
  84.                 c = a * b;
  85.                 c = c / a;
  86.                 c = a * b;
  87.                 c = c / a;
  88.                 c = a * b;
  89.                 c = c / a;
  90.                 }
  91. /*      printf ("Done\n"); */
  92.         }
  93.  
  94. #undef COUNT
  95.  
  96. /*------------------------------------------------------------------*/
  97.  
  98. /***************/
  99. /*    Flt.c    */
  100. /***************/
  101. flt() {
  102.    long i,j,loops; double a,b,c;
  103.    loops = 256000;
  104.  
  105.    for( i=1 ; i<=loops ; i++ ) {
  106.        j = loops - i;
  107.        a = (double)i;
  108.        b = (double)j;
  109.        c = b / a;
  110.        a = b - c;
  111.        b = c * a;
  112.        c = b + a;
  113.        }
  114.  
  115.    a = c + b;
  116.  
  117. }
  118.  
  119. /*------------------------------------------------------------------*/
  120.  
  121. /*
  122. ** savage.c -- floating point speed and accuracy test.  C version
  123. ** derived from BASIC version which appeared in Dr. Dobb's Journal,
  124. ** Sep. 1983, pp. 120-122.
  125. */
  126.  
  127. #define    ILOOP    25000
  128.  
  129. extern    double    tan(), atan(), exp(), log(), sqrt();
  130.  
  131. savage()
  132. {
  133. int i;
  134. double a;
  135.  
  136. /*    printf("start\n"); */
  137.     a = 1.0;
  138.     for (i = 1; i <= (ILOOP - 1); i++)
  139.         a = tan(atan(exp(log(sqrt(a*a))))) + 1.0;
  140. /*    printf("a = %20.14e\n", a); */
  141. /*    printf("done\n"); */
  142. }
  143.  
  144. /*------------------------------------------------------------------*/
  145.  
  146. /**************************************************************************
  147.  *                                                                        *
  148.  *      Whetstone benchmark in C.  This program is a translation of the   *
  149.  *      original Algol version in "A Synthetic Benchmark" by H.J. Curnow  *
  150.  *      and B.A. Wichman in Computer Journal, Vol  19 #1, February 1976.  *
  151.  *                                                                        *
  152.  *      Used to test compiler efficiency, optimization, and double        *
  153.  *      precision floating-point performance.  This version is specific   *
  154.  *      to the Turbo-Amiga and Amiga but it can be easily adapted to      *
  155.  *      other systems by replacing the clock() routine with your own. *
  156.  *                                                                        *
  157.  **************************************************************************/
  158.  
  159. #define ITERATIONS   10       /* 1 Million Whetstone instructions    */
  160.  
  161. /* #define POUT  */           /* Leave as is for 'Official' result.  */
  162.  
  163. /* #define MTEST */           /* define for Module timing results    */
  164.                               /* only.  Leave as is for 'Official'   */
  165.                               /* result.                             */
  166.  
  167. static double   D_e1[4];
  168. static double   D_t, D_t1, D_t2;
  169. static long     j, k, l;
  170.  
  171. whetd() {
  172.  
  173.    long        i;
  174.    long     n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11;
  175.    long     m, loops;
  176.    double   D_x, D_y, D_z, D_x1, D_x2, D_x3, D_x4;
  177.    
  178.    /************************/
  179.    /* initialize constants */
  180.    /************************/
  181.  
  182.    D_t   =   0.499975;
  183.    D_t1  =   0.500250;
  184.    D_t2  =   2.0;
  185.  
  186.    /***********************/
  187.    /* Set Module Weights. */
  188.    /***********************/
  189.  
  190.    m = 10;                    /* m = 10 is used to obtain better timing  */
  191.    loops = m * ITERATIONS;    /* accuracy only.  Slow systems should use */
  192.    n1  =   0 * loops;         /* m = 1.                                  */
  193.    n2  =  12 * loops;
  194.    n3  =  14 * loops;
  195.    n4  = 345 * loops;
  196.    n5  =   0 * loops;
  197.    n6  = 210 * loops;
  198.    n7  =  32 * loops;
  199.    n8  = 899 * loops;
  200.    n9  = 616 * loops;
  201.    n10 =   0 * loops;
  202.    n11 =  93 * loops;
  203.  
  204.    /*********************************/
  205.    /* MODULE 1:  simple identifiers */
  206.    /*********************************/
  207.  
  208.    D_x1 =  1.0;
  209.    D_x2 = -1.0;
  210.    D_x3 = -1.0;
  211.    D_x4 = -1.0;
  212.  
  213.    if( n1 > 0 )
  214.    {
  215.      for(i = 1; i <= n1; i++)
  216.      {
  217.       D_x1 = ( D_x1 + D_x2 + D_x3 - D_x4 ) * D_t;
  218.       D_x2 = ( D_x1 + D_x2 - D_x3 - D_x4 ) * D_t;
  219.       D_x3 = ( D_x1 - D_x2 + D_x3 + D_x4 ) * D_t;
  220.       D_x4 = (-D_x1 + D_x2 + D_x3 + D_x4 ) * D_t;
  221.      }
  222.    }
  223.  
  224.    /*****************************/
  225.    /* MODULE 2:  Array Elements */
  226.    /*****************************/
  227.    D_e1[0] =  1.0;        /* Start at element 0 in C, vice 1 in Fortran */
  228.    D_e1[1] = -1.0;
  229.    D_e1[2] = -1.0;
  230.    D_e1[3] = -1.0;
  231.  
  232.    if( n2 > 0 )
  233.    {
  234.      for (i = 1; i <= n2; i++)
  235.      {
  236.       D_e1[0] = ( D_e1[0] + D_e1[1] + D_e1[2] - D_e1[3] ) * D_t;
  237.       D_e1[1] = ( D_e1[0] + D_e1[1] - D_e1[2] + D_e1[3] ) * D_t;
  238.       D_e1[2] = ( D_e1[0] - D_e1[1] + D_e1[2] + D_e1[3] ) * D_t;
  239.       D_e1[3] = (-D_e1[0] + D_e1[1] + D_e1[2] + D_e1[3] ) * D_t;
  240.      }
  241.    }
  242.  
  243.    /*********************************/
  244.    /* MODULE 3:  Array as Parameter */
  245.    /*********************************/
  246.    if( n3 > 0 )
  247.    {
  248.      for (i = 1; i <= n3; i++)
  249.      {
  250.       D_pa(D_e1);
  251.      }
  252.    }
  253.  
  254.    /********************************/
  255.    /* MODULE 4:  Conditional Jumps */
  256.    /********************************/
  257.    j = 1;
  258.  
  259.    if( n4 > 0 )
  260.    {
  261.      for (i = 1; i <= n4; i++)
  262.      {
  263.       if (j == 1)
  264.          j = 2;
  265.       else
  266.          j = 3;
  267.  
  268.       if (j > 2)
  269.          j = 0;
  270.       else
  271.          j = 1;
  272.  
  273.       if (j < 1 )
  274.          j = 1;
  275.       else
  276.          j = 0;
  277.      }
  278.    }
  279.  
  280.    /**********************/
  281.    /* MODULE 5:  Omitted */
  282.    /**********************/
  283.  
  284.  
  285.    /*********************************/
  286.    /* MODULE 6:  Integer Arithmetic */
  287.    /*********************************/
  288.    j = 1;
  289.    k = 2;
  290.    l = 3;
  291.  
  292.    if( n6 > 0 )
  293.    {
  294.      for (i = 1; i <= n6; i++)
  295.      {
  296.       j = j * (k - j) * (l -k);
  297.       k = l * k - (l - j) * k;
  298.       l = (l - k) * (k + j);
  299.  
  300.       D_e1[l - 2] = j + k + l;          /* Remember we started at D_e1[0]. */
  301.       D_e1[k - 2] = j * k * l;          /* l-2 in C, vice l-1 in Fortran */
  302.      }
  303.    }
  304.  
  305.    /**************************************/
  306.    /* MODULE 7:  Trigonometric Functions */
  307.    /**************************************/
  308.    D_x = 0.5;
  309.    D_y = 0.5;
  310.  
  311.    if( n7 > 0 )
  312.    {
  313.      for(i = 1; i <= n7; i++)
  314.      {
  315.       D_x = D_t * atan(D_t2*sin(D_x)*cos(D_x)/(cos(D_x+D_y)+cos(D_x-D_y)-1.0));
  316.       D_y = D_t * atan(D_t2*sin(D_y)*cos(D_y)/(cos(D_x+D_y)+cos(D_x-D_y)-1.0));
  317.      }
  318.    }
  319.  
  320.    /******************************/
  321.    /* MODULE 8:  Procedure Calls */
  322.    /******************************/
  323.  
  324.    D_x = 1.0;
  325.    D_y = 1.0;
  326.    D_z = 1.0;
  327.  
  328.    if( n8 > 0 )
  329.    {
  330.      for (i = 1; i <= n8; i++)
  331.      {
  332.       D_p3(D_x, D_y, &D_z);
  333.      }
  334.    }
  335.  
  336.    /*******************************/
  337.    /* MODULE 9:  Array References */
  338.    /*******************************/
  339.  
  340.    j = 1;
  341.    k = 2;
  342.    l = 3;
  343.  
  344.    D_e1[0] = 1.0;
  345.    D_e1[1] = 2.0;
  346.    D_e1[2] = 3.0;
  347.  
  348.    if( n9 > 0 )
  349.    {
  350.      for(i = 1; i <= n9; i++)
  351.      {
  352.       D_p0();
  353.      }
  354.    }
  355.  
  356.    /**********************************/
  357.    /* MODULE 10:  Integer Arithmetic */
  358.    /**********************************/
  359.  
  360.    j = 2;
  361.    k = 3;
  362.  
  363.    if( n10 > 0 )
  364.    {
  365.      for(i = 1; i <= n10; i++)
  366.      {
  367.       j = j + k;
  368.       k = j + k;
  369.       j = k - j;
  370.       k = k - j - j;
  371.      }
  372.    }
  373.  
  374.    /**********************************/
  375.    /* MODULE 11:  Standard Functions */
  376.    /**********************************/
  377.  
  378.    D_x = 0.75;
  379.  
  380.    if( n11 > 0 )
  381.    {
  382.      for(i = 1; i <= n11; i++)
  383.      {
  384.       D_x = sqrt( exp( log(D_x) / D_t1) );
  385.      }
  386.    }
  387.  
  388.    /**************************/
  389.    /* End of Whetstone Tests */
  390.    /**************************/
  391.  
  392. #if 0    /* Done elsewhere now. */
  393.    stoptime  = clock();
  394.    benchtime = stoptime - starttime - nulltime;
  395.    D_x1 = (double)benchtime/100.0;
  396.    printf("   Benchtime(sec) = %lf\n",D_x1);
  397.    D_x2 = 100.0 * (double)loops / D_x1;
  398.    KWhets = (long)D_x2;
  399.    printf("   KWhets/sec     = %ld\n\n",KWhets);
  400. #endif   
  401.  
  402. }
  403.  
  404. /*******************/
  405. /* Subroutine pa() */
  406. /*******************/
  407.  
  408. D_pa(e)                /* Exactly as in the Algol 60 version, but we */
  409.                      /* could do away with that 'goto'.            */
  410. double e[4];
  411.  
  412. {
  413.    int j;
  414.  
  415.    j = 0;
  416.      lab:
  417.    e[0] = (  e[0] + e[1] + e[2] - e[3] ) * D_t;
  418.    e[1] = (  e[0] + e[1] - e[2] + e[3] ) * D_t;
  419.    e[2] = (  e[0] - e[1] + e[2] + e[3] ) * D_t;
  420.    e[3] = ( -e[0] + e[1] + e[2] + e[3] ) / D_t2;
  421.    j ++;
  422.  
  423.    if (j < 6)
  424.       goto lab;
  425. }
  426.  
  427. /************************/
  428. /* Subroutine p3(x,y,z) */
  429. /************************/
  430.  
  431. D_p3(x, y, z)
  432.  
  433. double x, y, *z;
  434.  
  435. {
  436.    x  = D_t * (x + y);
  437.    y  = D_t * (x + y);
  438.    *z = (x + y) /D_t2;
  439. }
  440.  
  441. /*******************/
  442. /* Subroutine p0() */
  443. /*******************/
  444.  
  445. D_p0()
  446. {
  447.    D_e1[j] = D_e1[k];
  448.    D_e1[k] = D_e1[l];
  449.    D_e1[l] = D_e1[j];
  450. }
  451.  
  452. static float    S_e1[4];
  453. static float    S_t, S_t1, S_t2; 
  454.  
  455. /* A true single-precision whestone can only be achieved by using
  456.    ANSI prototypes.  Without them, doubles get passed to p3 and converted
  457.    to float upon entry to p3.
  458.    If your compiler doesn't support prototypes, you could instead
  459.    pass pointers to x and y, for a slight decrease in efficiency.
  460. */   
  461.  
  462. #ifndef NO_PROTOTYPES
  463. S_p3(float,float,float*);
  464. #endif
  465.  
  466. whets() {
  467.  
  468.    long        i;
  469.    long     n1, n2, n3, n4, n5, n6, n7, n8, n9, n10, n11;
  470.    long     m, loops;
  471.    float    S_x, S_y, S_z, S_x1, S_x2, S_x3, S_x4;
  472.    
  473.    /************************/
  474.    /* initialize constants */
  475.    /************************/
  476.  
  477.    S_t   =   0.499975;
  478.    S_t1  =   0.500250;
  479.    S_t2  =   2.0;
  480.  
  481.    /***********************/
  482.    /* Set Module Weights. */
  483.    /***********************/
  484.  
  485.    m = 10;                    /* m = 10 is used to obtain better timing  */
  486.    loops = m * ITERATIONS;    /* accuracy only.  Slow systems should use */
  487.    n1  =   0 * loops;         /* m = 1.                                  */
  488.    n2  =  12 * loops;
  489.    n3  =  14 * loops;
  490.    n4  = 345 * loops;
  491.    n5  =   0 * loops;
  492.    n6  = 210 * loops;
  493.    n7  =  32 * loops;
  494.    n8  = 899 * loops;
  495.    n9  = 616 * loops;
  496.    n10 =   0 * loops;
  497.    n11 =  93 * loops;
  498.  
  499.    /*********************************/
  500.    /* MODULE 1:  simple identifiers */
  501.    /*********************************/
  502.  
  503.    S_x1 =  1.0;
  504.    S_x2 = -1.0;
  505.    S_x3 = -1.0;
  506.    S_x4 = -1.0;
  507.  
  508.    if( n1 > 0 )
  509.    {
  510.      for(i = 1; i <= n1; i++)
  511.      {
  512.       S_x1 = ( S_x1 + S_x2 + S_x3 - S_x4 ) * S_t;
  513.       S_x2 = ( S_x1 + S_x2 - S_x3 - S_x4 ) * S_t;
  514.       S_x3 = ( S_x1 - S_x2 + S_x3 + S_x4 ) * S_t;
  515.       S_x4 = (-S_x1 + S_x2 + S_x3 + S_x4 ) * S_t;
  516.      }
  517.    }
  518.  
  519.    /*****************************/
  520.    /* MODULE 2:  Array Elements */
  521.    /*****************************/
  522.    S_e1[0] =  1.0;        /* Start at element 0 in C, vice 1 in Fortran */
  523.    S_e1[1] = -1.0;
  524.    S_e1[2] = -1.0;
  525.    S_e1[3] = -1.0;
  526.  
  527.    if( n2 > 0 )
  528.    {
  529.      for (i = 1; i <= n2; i++)
  530.      {
  531.       S_e1[0] = ( S_e1[0] + S_e1[1] + S_e1[2] - S_e1[3] ) * S_t;
  532.       S_e1[1] = ( S_e1[0] + S_e1[1] - S_e1[2] + S_e1[3] ) * S_t;
  533.       S_e1[2] = ( S_e1[0] - S_e1[1] + S_e1[2] + S_e1[3] ) * S_t;
  534.       S_e1[3] = (-S_e1[0] + S_e1[1] + S_e1[2] + S_e1[3] ) * S_t;
  535.      }
  536.    }
  537.  
  538.    /*********************************/
  539.    /* MODULE 3:  Array as Parameter */
  540.    /*********************************/
  541.    if( n3 > 0 )
  542.    {
  543.      for (i = 1; i <= n3; i++)
  544.      {
  545.       S_pa(S_e1);
  546.      }
  547.    }
  548.  
  549.    /********************************/
  550.    /* MODULE 4:  Conditional Jumps */
  551.    /********************************/
  552.    j = 1;
  553.  
  554.    if( n4 > 0 )
  555.    {
  556.      for (i = 1; i <= n4; i++)
  557.      {
  558.       if (j == 1)
  559.          j = 2;
  560.       else
  561.          j = 3;
  562.  
  563.       if (j > 2)
  564.          j = 0;
  565.       else
  566.          j = 1;
  567.  
  568.       if (j < 1 )
  569.          j = 1;
  570.       else
  571.          j = 0;
  572.      }
  573.    }
  574.  
  575.    /**********************/
  576.    /* MODULE 5:  Omitted */
  577.    /**********************/
  578.  
  579.  
  580.    /*********************************/
  581.    /* MODULE 6:  Integer Arithmetic */
  582.    /*********************************/
  583.    j = 1;
  584.    k = 2;
  585.    l = 3;
  586.  
  587.    if( n6 > 0 )
  588.    {
  589.      for (i = 1; i <= n6; i++)
  590.      {
  591.       j = j * (k - j) * (l -k);
  592.       k = l * k - (l - j) * k;
  593.       l = (l - k) * (k + j);
  594.  
  595.       S_e1[l - 2] = j + k + l;          /* Remember we started at S_e1[0]. */
  596.       S_e1[k - 2] = j * k * l;          /* l-2 in C, vice l-1 in Fortran */
  597.      }
  598.    }
  599.  
  600.    /**************************************/
  601.    /* MODULE 7:  Trigonometric Functions */
  602.    /**************************************/
  603.    S_x = 0.5;
  604.    S_y = 0.5;
  605.  
  606.    if( n7 > 0 )
  607.    {
  608.      for(i = 1; i <= n7; i++)
  609.      {
  610.       S_x = S_t * atan(S_t2*sin(S_x)*cos(S_x)/(cos(S_x+S_y)+cos(S_x-S_y)-1));
  611.       S_y = S_t * atan(S_t2*sin(S_y)*cos(S_y)/(cos(S_x+S_y)+cos(S_x-S_y)-1));
  612.      }
  613.    }
  614.  
  615.    /******************************/
  616.    /* MODULE 8:  Procedure Calls */
  617.    /******************************/
  618.  
  619.    S_x = 1.0;
  620.    S_y = 1.0;
  621.    S_z = 1.0;
  622.  
  623.    if( n8 > 0 )
  624.    {
  625.      for (i = 1; i <= n8; i++)
  626.      {
  627.       S_p3(S_x, S_y, &S_z);
  628.      }
  629.    }
  630.  
  631.    /*******************************/
  632.    /* MODULE 9:  Array References */
  633.    /*******************************/
  634.  
  635.    j = 1;
  636.    k = 2;
  637.    l = 3;
  638.  
  639.    S_e1[0] = 1.0;
  640.    S_e1[1] = 2.0;
  641.    S_e1[2] = 3.0;
  642.  
  643.    if( n9 > 0 )
  644.    {
  645.      for(i = 1; i <= n9; i++)
  646.      {
  647.       S_p0();
  648.      }
  649.    }
  650.  
  651.    /**********************************/
  652.    /* MODULE 10:  Integer Arithmetic */
  653.    /**********************************/
  654.  
  655.    j = 2;
  656.    k = 3;
  657.  
  658.    if( n10 > 0 )
  659.    {
  660.      for(i = 1; i <= n10; i++)
  661.      {
  662.       j = j + k;
  663.       k = j + k;
  664.       j = k - j;
  665.       k = k - j - j;
  666.      }
  667.    }
  668.  
  669.    /**********************************/
  670.    /* MODULE 11:  Standard Functions */
  671.    /**********************************/
  672.  
  673.    S_x = 0.75;
  674.  
  675.    if( n11 > 0 )
  676.    {
  677.      for(i = 1; i <= n11; i++)
  678.      {
  679.       S_x = sqrt( exp( log(S_x) / S_t1) );
  680.      }
  681.    }
  682.  
  683.    /**************************/
  684.    /* End of Whetstone Tests */
  685.    /**************************/
  686.  
  687. #if 0    /* Done elsewhere now. */
  688.    stoptime  = clock();
  689.    benchtime = stoptime - starttime - nulltime;
  690.    S_x1 = (double)benchtime/100.0;
  691.    printf("   Benchtime(sec) = %lf\n",S_x1);
  692.    S_x2 = 100.0 * (double)loops / S_x1;
  693.    KWhets = (long)S_x2;
  694.    printf("   KWhets/sec     = %ld\n\n",KWhets);
  695. #endif   
  696.  
  697. }
  698.  
  699. /*******************/
  700. /* Subroutine pa() */
  701. /*******************/
  702.  
  703. S_pa(e)                /* Exactly as in the Algol 60 version, but we */
  704.                      /* could do away with that 'goto'.            */
  705. float e[4];
  706.  
  707. {
  708.    int j;
  709.  
  710.    j = 0;
  711.      lab:
  712.    e[0] = (  e[0] + e[1] + e[2] - e[3] ) * S_t;
  713.    e[1] = (  e[0] + e[1] - e[2] + e[3] ) * S_t;
  714.    e[2] = (  e[0] - e[1] + e[2] + e[3] ) * S_t;
  715.    e[3] = ( -e[0] + e[1] + e[2] + e[3] ) / S_t2;
  716.    j ++;
  717.  
  718.    if (j < 6)
  719.       goto lab;
  720. }
  721.  
  722. /************************/
  723. /* Subroutine p3(x,y,z) */
  724. /************************/
  725.  
  726. S_p3(x, y, z)
  727.  
  728. float x, y, *z;
  729.  
  730. {
  731.    x  = S_t * (x + y);
  732.    y  = S_t * (x + y);
  733.    *z = (x + y) /S_t2;
  734. }
  735.  
  736. /*******************/
  737. /* Subroutine p0() */
  738. /*******************/
  739.  
  740. S_p0()
  741. {
  742.    S_e1[j] = S_e1[k];
  743.    S_e1[k] = S_e1[l];
  744.    S_e1[l] = S_e1[j];
  745. }
  746.  
  747. /*-- End Of Whetstone C Source Code -----------------*/
  748. /*------------------------------------------------------------------*/
  749. /* Autodesk benchmark.                            */
  750. /*------------------------------------------------------------------*/
  751.  
  752. /*
  753.  
  754.     This is a complete optical design raytracing algorithm,
  755.     stripped of its user interface and recast into portable C.  It
  756.     not only determines execution speed on an extremely floating
  757.     point (including trig function) intensive real-world
  758.     application, it checks accuracy on an algorithm that is
  759.     exquisitely sensitive to errors.  The performance of this
  760.     program is typically far more sensitive to changes in the
  761.     efficiency of the trigonometric library routines than the
  762.     average floating point program.
  763.  
  764.     The benchmark may be compiled in two modes.  If the symbol
  765.     INTRIG is defined, built-in trigonometric and square root
  766.     routines will be used for all calculations.  Timings made with
  767.         INTRIG defined reflect the machine's basic floating point
  768.     performance for the arithmetic operators.  If INTRIG is not
  769.     defined, the system library <math.h> functions are used.  Results
  770.         with INTRIG not defined reflect the system's library performance
  771.     and/or floating point hardware support for trig functions and
  772.     square root.  Results with INTRIG defined are a good guide to
  773.     general floating point performance, while results with INTRIG
  774.     undefined indicate the performance of an application which is
  775.     math function intensive.
  776.  
  777.     Special note regarding errors in accuracy: this program has
  778.     generated numbers identical to the last digit it formats and
  779.     checks on the following machines, floating point architectures,
  780.     and languages:
  781.  
  782.     IBM PC / XT / AT  Lattice C IEEE 64 bit, 80 bit temporaries
  783.               High C    same, in line 80x87 code
  784.                           BASICA    "Double precision"
  785.               Quick BASIC IEEE double precision, software routines
  786.  
  787.     Sun 3          C        IEEE 64 bit, 80 bit temporaries,
  788.                     in-line 68881 code, in-line FPA code.
  789.  
  790.         MicroVAX II       C         Vax "G" format floating point
  791.  
  792.     Macintosh Plus      MPW C     SANE floating point, IEEE 64 bit format
  793.                     implemented in ROM.
  794.  
  795.     Inaccuracies reported by this program should be taken VERY
  796.     SERIOUSLY INDEED, as the program has been demonstrated to be
  797.     invariant under changes in floating point format, as long as
  798.     the format is a recognised double precision format.  If you
  799.     encounter errors, please remember that they are just as likely
  800.     to be in the floating point editing library or the
  801.     trigonometric libraries as in the low level operator code.
  802.  
  803.     The benchmark assumes that results are basically reliable, and
  804.     only tests the last result computed against the reference.  If
  805.         you're running on a suspect system you can compile this program
  806.     with ACCURACY defined.    This will generate a version which
  807.     executes as an infinite loop, performing the ray trace and
  808.     checking the results on every pass.  All incorrect results will
  809.     be reported.
  810.  
  811.     Representative timings are given below.  All have been normalised
  812.     as if run for 1000 iterations.
  813.  
  814.   Time in seconds           Computer, Compiler, and notes
  815.  Normal      INTRIG
  816.  
  817.  3466.00    4031.00    Commodore 128, 2 Mhz 8510 with software floating
  818.             point.    Abacus Software/Data-Becker Super-C 128,
  819.             version 3.00, run in fast (2 Mhz) mode.  Note:
  820.             the results generated by this system differed
  821.             from the reference results in the 8th to 10th
  822.             decimal place.
  823.  
  824.  3290.00        IBM PC/AT 6 Mhz, Microsoft/IBM BASICA version A3.00.
  825.                         Run with the "/d" switch, software floating point.
  826.  
  827.  2131.50        IBM PC/AT 6 Mhz, Lattice C version 2.14, small model.
  828.             This version of Lattice compiles subroutine
  829.             calls which either do software floating point
  830.             or use the 80x87.  The machine on which I ran
  831.             this had an 80287, but the results were so bad
  832.             I wonder if it was being used.
  833.  
  834.  1598.00        Macintosh Plus, MPW C, SANE Software floating point.
  835.  
  836.   404.00        IBM PC/AT 6 Mhz, Microsoft QuickBASIC version 2.0.
  837.             Software floating point.
  838.  
  839.   165.15        IBM PC/AT 6 Mhz, Metaware High C version 1.3, small
  840.             model.    This was compiled to call subroutines for
  841.             floating point, and the machine contained an 80287
  842.             which was used by the subroutines.
  843.  
  844.   143.20        Macintosh II, MPW C, SANE calls.  I was unable to
  845.             determine whether SANE was using the 68881 chip or
  846.             not.
  847.  
  848.   121.80        Sun 3/160 16 Mhz, Sun C.  Compiled with -fsoft switch
  849.             which executes floating point in software.
  850.  
  851.    78.78        IBM RT PC (Model 6150).  IBM AIX 1.0 C compiler
  852.             with -O switch.
  853.  
  854.    66.36     206.35    IBM PC/AT 6Mhz, Metaware High C version 1.3, small
  855.             model.    This was compiled with in-line code for the
  856.             80287 math coprocessor.  Trig functions still call
  857.             library routines.
  858.  
  859.    17.18        Apollo DN-3000, 12 Mhz 68020 with 68881, compiled
  860.             with in-line code for the 68881 coprocessor.
  861.             According to Apollo, the library routines are chosen
  862.             at runtime based on coprocessor presence.  Since the
  863.             coprocessor was present, the library is supposed to
  864.             use in-line floating point code.
  865.  
  866.    15.55      27.56    VAXstation II GPX.  Compiled and executed under
  867.             VAX/VMS C.
  868.  
  869.    15.14      37.93    Macintosh II, Unix system V.  Green Hills 68020
  870.             Unix compiler with in-line code for the 68881
  871.             coprocessor (-O -ZI switches).
  872.  
  873.    12.69        Sun 3/160 16 Mhz, Sun C.  Compiled with -fswitch,
  874.             which calls a subroutine to select the fastest
  875.             floating point processor.  This was using the 68881.
  876.  
  877.    11.74      26.73    Compaq Deskpro 386, 16 Mhz 80386 with 16 Mhz 80387.
  878.             Metaware High C version 1.3, compiled with in-line
  879.             for the math coprocessor (but not optimised for the
  880.             80386/80387).  Trig functions still call library
  881.             routines.
  882.  
  883.     8.43      30.49    Sun 3/160 16 Mhz, Sun C.  Compiled with -f68881,
  884.             generating in-line MC68881 instructions.  Trig
  885.             functions still call library routines.
  886.  
  887.     6.29      25.17    Sun 3/260 25 Mhz, Sun C.  Compiled with -f68881,
  888.             generating in-line MC68881 instructions.  Trig
  889.             functions still call library routines.
  890.  
  891.     6.00      23.00    Compaq Deskpro 386, 16 Mhz 80386 with 16 Mhz 80387.
  892.             Metaware High C version 1.3, compiled with in-line
  893.             for the math coprocessor (optimised for the
  894.             80386/80387).   
  895.  
  896.     5.9       23.00    Compaq Deskpro 386, 16 Mhz 80386 with 16 Mhz 80387.
  897.             Metaware High C version 1.4, compiled with in-line
  898.             for the math coprocessor (optimised for the
  899.             80386/80387).    
  900.  
  901.     5.25      ---       Intel 386/20, 16 Mhz 80386 with 16 Mhz 80387.
  902.             Metaware High C version 1.4, compiled with in-line
  903.             for the math coprocessor (optimised for the
  904.             80386/80387).    
  905.  
  906. */
  907.  
  908. #define cot(x) (1.0 / tan(x))
  909. #define INTRIG_cot(x) (1.0 / INTRIG_tan(x))
  910.  
  911. #define max_surfaces 10
  912.  
  913. /*  Local variables  */
  914.  
  915. static char tbfr[132];
  916.  
  917. static short current_surfaces;
  918. static short paraxial;
  919.  
  920. static double clear_aperture;
  921.  
  922. static double aberr_lspher;
  923. static double aberr_osc;
  924. static double aberr_lchrom;
  925.  
  926. static double max_lspher;
  927. static double max_osc;
  928. static double max_lchrom;
  929.  
  930. static double radius_of_curvature;
  931. static double object_distance;
  932. static double ray_height;
  933. static double axis_slope_angle;
  934. static double from_index;
  935. static double to_index;
  936.  
  937. static double spectral_line[9];
  938. static double s[max_surfaces][5];
  939. static double od_sa[2][2];
  940.  
  941. static char outarr[8][80];       /* Computed output of program goes here */
  942.  
  943. int itercount;               /* The iteration counter for the main loop
  944.                       in the program is made global so that
  945.                       the compiler should not be allowed to
  946.                       optimise out the loop over the ray
  947.                       tracing code. */
  948.  
  949. /* The test case used in this program is the design for a 4 inch
  950.    achromatic telescope objective used as the example in Wyld's
  951.    classic work on ray tracing by hand, given in Amateur Telescope
  952.    Making, Volume 3. */
  953.  
  954. static double testcase[4][4] = {
  955.     {27.05, 1.5137, 63.6, 0.52},
  956.     {-16.68, 1, 0, 0.138},
  957.     {-16.68, 1.6164, 36.7, 0.38},
  958.     {-78.1, 1, 0, 0}
  959. };
  960.  
  961. /*  Internal trig functions (used only if INTRIG is defined).  These
  962.     standard functions may be enabled to obtain timings that reflect
  963.     the machine's floating point performance rather than the speed of
  964.     its trig function evaluation.  */
  965.  
  966. #define INTRIG_fabs(x)  ((x < 0.0) ? -x : x)
  967.  
  968. #define pic 3.1415926535897932
  969.  
  970. /*  Commonly used constants  */
  971.  
  972. static double pi = pic,
  973.     twopi =pic * 2.0,
  974.     piover4 = pic / 4.0,
  975.     fouroverpi = 4.0 / pic,
  976.     piover2 = pic / 2.0;
  977.  
  978. /*  Coefficients for ATAN evaluation  */
  979.  
  980. static double atanc[] = {
  981.     0.0,
  982.     0.4636476090008061165,
  983.     0.7853981633974483094,
  984.     0.98279372324732906714,
  985.     1.1071487177940905022,
  986.     1.1902899496825317322,
  987.     1.2490457723982544262,
  988.     1.2924966677897852673,
  989.     1.3258176636680324644
  990. };
  991.  
  992. /*  aint(x)      Return integer part of number.  Truncates towards 0     */
  993.  
  994. double aint(x)
  995. double x;
  996. {
  997.     long l;
  998.  
  999.     /*  Note that this routine cannot handle the full floating point
  1000.         number range.  This function should be in the machine-dependent
  1001.         floating point library!  */
  1002.  
  1003.     l = x;
  1004.     if ((int)(-0.5) != 0  &&  l < 0 )
  1005.        l++;
  1006.     x = l;
  1007.     return x;
  1008. }
  1009.  
  1010. /*  sin(x)      Return sine, x in radians  */
  1011.  
  1012. double INTRIG_sin(x)
  1013. double x;
  1014. {
  1015.     int sign;
  1016.     double y, r, z;
  1017.  
  1018.     x = (((sign= (x < 0.0)) != 0) ? -x: x);
  1019.  
  1020.     if (x > twopi)
  1021.        x -= (aint(x / twopi) * twopi);
  1022.  
  1023.     if (x > pi) {
  1024.        x -= pi;
  1025.        sign = !sign;
  1026.     }
  1027.  
  1028.     if (x > piover2)
  1029.        x = pi - x;
  1030.  
  1031.     if (x < piover4) {
  1032.        y = x * fouroverpi;
  1033.        z = y * y;
  1034.        r = y * (((((((-0.202253129293E-13 * z + 0.69481520350522E-11) * z -
  1035.           0.17572474176170806E-8) * z + 0.313361688917325348E-6) * z -
  1036.           0.365762041821464001E-4) * z + 0.249039457019271628E-2) * z -
  1037.           0.0807455121882807815) * z + 0.785398163397448310);
  1038.     } else {
  1039.        y = (piover2 - x) * fouroverpi;
  1040.        z = y * y;
  1041.        r = ((((((-0.38577620372E-12 * z + 0.11500497024263E-9) * z -
  1042.           0.2461136382637005E-7) * z + 0.359086044588581953E-5) * z -
  1043.           0.325991886926687550E-3) * z + 0.0158543442438154109) * z -
  1044.           0.308425137534042452) * z + 1.0;
  1045.     }
  1046.     return sign ? -r : r;
  1047. }
  1048.  
  1049. /*  cos(x)      Return cosine, x in radians, by identity  */
  1050.  
  1051. double INTRIG_cos(x)
  1052. double x;
  1053. {
  1054.     x = (x < 0.0) ? -x : x;
  1055.     if (x > twopi)              /* Do range reduction here to limit */
  1056.        x = x - (aint(x / twopi) * twopi); /* roundoff on add of PI/2    */
  1057.     return INTRIG_sin(x + piover2);
  1058. }
  1059.  
  1060. /*  tan(x)      Return tangent, x in radians, by identity  */
  1061.  
  1062. static double INTRIG_tan(x)
  1063. double x;
  1064. {
  1065.     return INTRIG_sin(x) / INTRIG_cos(x);
  1066. }
  1067.  
  1068. /*  sqrt(x)      Return square root.  Initial guess, then Newton-
  1069.           Raphson refinement  */
  1070.  
  1071. double INTRIG_sqrt(x)
  1072. double x;
  1073. {
  1074.     double c, cl, y;
  1075.     int n;
  1076.  
  1077.     if (x == 0.0)
  1078.        return 0.0;
  1079.  
  1080.     if (x < 0.0) {
  1081.        fprintf(stderr,
  1082.               "\nGood work!  You tried to take the square root of %g",
  1083.          x);
  1084.        fprintf(stderr,
  1085.               "\nunfortunately, that is too complex for me to handle.\n");
  1086.        exit(1);
  1087.     }
  1088.  
  1089.     y = (0.154116 + 1.893872 * x) / (1.0 + 1.047988 * x);
  1090.  
  1091.     c = (y - x / y) / 2.0;
  1092.     cl = 0.0;
  1093.     for (n = 50; c != cl && n--;) {
  1094.        y = y - c;
  1095.        cl = c;
  1096.        c = (y - x / y) / 2.0;
  1097.     }
  1098.     return y;
  1099. }
  1100.  
  1101. /*  atan(x)      Return arctangent in radians,
  1102.           range -pi/2 to pi/2  */
  1103.  
  1104. double INTRIG_atan(x)
  1105. double x;
  1106. {
  1107.     int sign, l, y;
  1108.     double a, b, z;
  1109.  
  1110.     x = (((sign = (x < 0.0)) != 0) ? -x : x);
  1111.     l = 0;
  1112.  
  1113.     if (x >= 4.0) {
  1114.        l = -1;
  1115.        x = 1.0 / x;
  1116.        y = 0;
  1117.        goto atl;
  1118.     } else {
  1119.        if (x < 0.25) {
  1120.           y = 0;
  1121.           goto atl;
  1122.        }
  1123.     }
  1124.  
  1125.     y = aint(x / 0.5);
  1126.     z = y * 0.5;
  1127.     x = (x - z) / (x * z + 1);
  1128.  
  1129. atl:
  1130.     z = x * x;
  1131.     b = ((((893025.0 * z + 49116375.0) * z + 425675250.0) * z +
  1132.         1277025750.0) * z + 1550674125.0) * z + 654729075.0;
  1133.     a = (((13852575.0 * z + 216602100.0) * z + 891080190.0) * z +
  1134.         1332431100.0) * z + 654729075.0;
  1135.     a = (a / b) * x + atanc[y];
  1136.     if (l)
  1137.        a=piover2 - a;
  1138.     return sign ? -a : a;
  1139. }
  1140.  
  1141. /*  atan2(y,x)      Return arctangent in radians of y/x,
  1142.           range -pi to pi  */
  1143.  
  1144. double INTRIG_atan2(y, x)
  1145. double y, x;
  1146. {
  1147.     double temp;
  1148.  
  1149.     if (x == 0.0) {
  1150.        if (y == 0.0)   /*  Special case: atan2(0,0) = 0  */
  1151.           return 0.0;
  1152.        else if (y > 0)
  1153.           return piover2;
  1154.        else
  1155.           return -piover2;
  1156.     }
  1157.     temp = INTRIG_atan(y / x);
  1158.     if (x < 0.0) {
  1159.        if (y >= 0.0)
  1160.           temp += pic;
  1161.        else
  1162.           temp -= pic;
  1163.     }
  1164.     return temp;
  1165. }
  1166.  
  1167. /*  asin(x)      Return arcsine in radians of x  */
  1168.  
  1169. double INTRIG_asin(x)
  1170. double x;
  1171. {
  1172.     if (INTRIG_fabs(x)>1.0) {
  1173.        fprintf(stderr,
  1174.               "\nInverse trig functions lose much of their gloss when");
  1175.        fprintf(stderr,
  1176.               "\ntheir arguments are greater than 1, such as the");
  1177.        fprintf(stderr,
  1178.               "\nvalue %g you passed.\n", x);
  1179.        exit(1);
  1180.     }
  1181.     return INTRIG_atan2(x, INTRIG_sqrt(1 - x * x));
  1182. }
  1183.  
  1184. /*  Perform ray trace in specific spectral line  */
  1185.  
  1186. static trace_line(line, ray_h)
  1187. int line;
  1188. double ray_h;
  1189. {
  1190.     int i;
  1191.  
  1192.     object_distance = 0.0;
  1193.     ray_height = ray_h;
  1194.     from_index = 1.0;
  1195.  
  1196.     for (i = 1; i <= current_surfaces; i++) {
  1197.        radius_of_curvature = s[i][1];
  1198.        to_index = s[i][2];
  1199.        if (to_index > 1.0)
  1200.           to_index = to_index + ((spectral_line[4] -
  1201.          spectral_line[line]) /
  1202.          (spectral_line[3] - spectral_line[6])) * ((s[i][2] - 1.0) /
  1203.          s[i][3]);
  1204.        transit_surface();
  1205.        from_index = to_index;
  1206.        if (i < current_surfaces)
  1207.           object_distance = object_distance - s[i][4];
  1208.     }
  1209. }
  1210.  
  1211. static INTRIG_trace_line(line, ray_h)
  1212. int line;
  1213. double ray_h;
  1214. {
  1215.     int i;
  1216.  
  1217.     object_distance = 0.0;
  1218.     ray_height = ray_h;
  1219.     from_index = 1.0;
  1220.  
  1221.     for (i = 1; i <= current_surfaces; i++) {
  1222.        radius_of_curvature = s[i][1];
  1223.        to_index = s[i][2];
  1224.        if (to_index > 1.0)
  1225.           to_index = to_index + ((spectral_line[4] -
  1226.          spectral_line[line]) /
  1227.          (spectral_line[3] - spectral_line[6])) * ((s[i][2] - 1.0) /
  1228.          s[i][3]);
  1229.        INTRIG_transit_surface();
  1230.        from_index = to_index;
  1231.        if (i < current_surfaces)
  1232.           object_distance = object_distance - s[i][4];
  1233.     }
  1234. }
  1235.  
  1236. static transit_surface() {
  1237.     double iang,           /* Incidence angle */
  1238.            rang,           /* Refraction angle */
  1239.            iang_sin,       /* Incidence angle sin */
  1240.            rang_sin,       /* Refraction angle sin */
  1241.            old_axis_slope_angle, sagitta;
  1242.  
  1243.     if (paraxial) {
  1244.        if (radius_of_curvature != 0.0) {
  1245.           if (object_distance == 0.0) {
  1246.          axis_slope_angle = 0.0;
  1247.          iang_sin = ray_height / radius_of_curvature;
  1248.           } else
  1249.          iang_sin = ((object_distance -
  1250.             radius_of_curvature) / radius_of_curvature) *
  1251.             axis_slope_angle;
  1252.  
  1253.           rang_sin = (from_index / to_index) *
  1254.          iang_sin;
  1255.           old_axis_slope_angle = axis_slope_angle;
  1256.           axis_slope_angle = axis_slope_angle +
  1257.          iang_sin - rang_sin;
  1258.           if (object_distance != 0.0)
  1259.          ray_height = object_distance * old_axis_slope_angle;
  1260.           object_distance = ray_height / axis_slope_angle;
  1261.           return;
  1262.        }
  1263.        object_distance = object_distance * (to_index / from_index);
  1264.        axis_slope_angle = axis_slope_angle * (from_index / to_index);
  1265.        return;
  1266.     }
  1267.  
  1268.     if (radius_of_curvature != 0.0) {
  1269.        if (object_distance == 0.0) {
  1270.           axis_slope_angle = 0.0;
  1271.           iang_sin = ray_height / radius_of_curvature;
  1272.        } else {
  1273.           iang_sin = ((object_distance -
  1274.          radius_of_curvature) / radius_of_curvature) *
  1275.          sin(axis_slope_angle);
  1276.        }
  1277.        iang = asin(iang_sin);
  1278.        rang_sin = (from_index / to_index) * iang_sin;
  1279.        old_axis_slope_angle = axis_slope_angle;
  1280.        axis_slope_angle = axis_slope_angle + iang - asin(rang_sin);
  1281.        sagitta = sin((old_axis_slope_angle + iang) / 2.0);
  1282.        sagitta = 2.0 * radius_of_curvature*sagitta*sagitta;
  1283.        object_distance = ((radius_of_curvature * sin(
  1284.           old_axis_slope_angle + iang)) *
  1285.           cot(axis_slope_angle)) + sagitta;
  1286.        return;
  1287.     }
  1288.  
  1289.     rang = -asin((from_index / to_index) *
  1290.        sin(axis_slope_angle));
  1291.     object_distance = object_distance * ((to_index *
  1292.        cos(-rang)) / (from_index *
  1293.        cos(axis_slope_angle)));
  1294.     axis_slope_angle = -rang;
  1295. }
  1296.  
  1297. static INTRIG_transit_surface() {
  1298.     double iang,           /* Incidence angle */
  1299.            rang,           /* Refraction angle */
  1300.            iang_sin,       /* Incidence angle sin */
  1301.            rang_sin,       /* Refraction angle sin */
  1302.            old_axis_slope_angle, sagitta;
  1303.  
  1304.     if (paraxial) {
  1305.        if (radius_of_curvature != 0.0) {
  1306.           if (object_distance == 0.0) {
  1307.          axis_slope_angle = 0.0;
  1308.          iang_sin = ray_height / radius_of_curvature;
  1309.           } else
  1310.          iang_sin = ((object_distance -
  1311.             radius_of_curvature) / radius_of_curvature) *
  1312.             axis_slope_angle;
  1313.  
  1314.           rang_sin = (from_index / to_index) * iang_sin;
  1315.           old_axis_slope_angle = axis_slope_angle;
  1316.           axis_slope_angle = axis_slope_angle + iang_sin - rang_sin;
  1317.           if (object_distance != 0.0)
  1318.          ray_height = object_distance * old_axis_slope_angle;
  1319.           object_distance = ray_height / axis_slope_angle;
  1320.           return;
  1321.        }
  1322.        object_distance = object_distance * (to_index / from_index);
  1323.        axis_slope_angle = axis_slope_angle * (from_index / to_index);
  1324.        return;
  1325.     }
  1326.  
  1327.     if (radius_of_curvature != 0.0) {
  1328.        if (object_distance == 0.0) {
  1329.           axis_slope_angle = 0.0;
  1330.           iang_sin = ray_height / radius_of_curvature;
  1331.        } else {
  1332.           iang_sin = ((object_distance -
  1333.          radius_of_curvature) / radius_of_curvature) *
  1334.          INTRIG_sin(axis_slope_angle);
  1335.        }
  1336.        iang = INTRIG_asin(iang_sin);
  1337.        rang_sin = (from_index / to_index) * iang_sin;
  1338.        old_axis_slope_angle = axis_slope_angle;
  1339.        axis_slope_angle = axis_slope_angle +
  1340.           iang - INTRIG_asin(rang_sin);
  1341.        sagitta = INTRIG_sin((old_axis_slope_angle + iang) / 2.0);
  1342.        sagitta = 2.0 * radius_of_curvature*sagitta*sagitta;
  1343.        object_distance = ((radius_of_curvature * INTRIG_sin(
  1344.           old_axis_slope_angle + iang)) *
  1345.           INTRIG_cot(axis_slope_angle)) + sagitta;
  1346.        return;
  1347.     }
  1348.  
  1349.     rang = -INTRIG_asin((from_index / to_index) *
  1350.        INTRIG_sin(axis_slope_angle));
  1351.     object_distance = object_distance * ((to_index *
  1352.        INTRIG_cos(-rang)) / (from_index *
  1353.        INTRIG_cos(axis_slope_angle)));
  1354.     axis_slope_angle = -rang;
  1355. }
  1356.  
  1357. /*  Initialise when called the first time  */
  1358.  
  1359. autodesk_init() {
  1360.     int i, j;
  1361.     spectral_line[1] = 7621.0;     /* A */
  1362.     spectral_line[2] = 6869.955;     /* B */
  1363.     spectral_line[3] = 6562.816;     /* C */
  1364.     spectral_line[4] = 5895.944;     /* D */
  1365.     spectral_line[5] = 5269.557;     /* E */
  1366.     spectral_line[6] = 4861.344;     /* F */
  1367.         spectral_line[7] = 4340.477;     /* G'*/
  1368.     spectral_line[8] = 3968.494;     /* H */
  1369.  
  1370.     /* Load test case into working array */
  1371.  
  1372.     clear_aperture = 4.0;
  1373.     current_surfaces = 4;
  1374.     for (i = 0; i < current_surfaces; i++)
  1375.        for (j = 0; j < 4; j++)
  1376.           s[i + 1][j + 1] = testcase[i][j];
  1377.  
  1378.        }
  1379.  
  1380. autodesk() {
  1381.     double od_fline, od_cline;
  1382.     int niter = 1000;           /* Iteration counter */
  1383.     
  1384.     /* Perform ray trace the specified number of times. */
  1385.  
  1386.     for (itercount = 0; itercount < niter; itercount++) {
  1387.  
  1388.        for (paraxial = 0; paraxial <= 1; paraxial++) {
  1389.  
  1390.           /* Do main trace in D light */
  1391.  
  1392.           trace_line(4, clear_aperture / 2.0);
  1393.           od_sa[paraxial][0] = object_distance;
  1394.           od_sa[paraxial][1] = axis_slope_angle;
  1395.        }
  1396.        paraxial = FALSE;
  1397.  
  1398.        /* Trace marginal ray in C */
  1399.  
  1400.        trace_line(3, clear_aperture / 2.0);
  1401.        od_cline = object_distance;
  1402.  
  1403.        /* Trace marginal ray in F */
  1404.  
  1405.        trace_line(6, clear_aperture / 2.0);
  1406.        od_fline = object_distance;
  1407.  
  1408.        aberr_lspher = od_sa[1][0] - od_sa[0][0];
  1409.        aberr_osc = 1.0 - (od_sa[1][0] * od_sa[1][1]) /
  1410.           (sin(od_sa[0][1]) * od_sa[0][0]);
  1411.        aberr_lchrom = od_fline - od_cline;
  1412.        max_lspher = sin(od_sa[0][1]);
  1413.  
  1414.        /* D light */
  1415.  
  1416.        max_lspher = 0.0000926 / (max_lspher * max_lspher);
  1417.        max_osc = 0.0025;
  1418.        max_lchrom = max_lspher;
  1419.     }
  1420.  
  1421. }
  1422.  
  1423. INTRIG_autodesk() {
  1424.     double od_fline, od_cline;
  1425.     int niter = 1000;           /* Iteration counter */
  1426.     
  1427.     /* Perform ray trace the specified number of times. */
  1428.  
  1429.     for (itercount = 0; itercount < niter; itercount++) {
  1430.  
  1431.        for (paraxial = 0; paraxial <= 1; paraxial++) {
  1432.  
  1433.           /* Do main trace in D light */
  1434.  
  1435.           INTRIG_trace_line(4, clear_aperture / 2.0);
  1436.           od_sa[paraxial][0] = object_distance;
  1437.           od_sa[paraxial][1] = axis_slope_angle;
  1438.        }
  1439.        paraxial = FALSE;
  1440.  
  1441.        /* Trace marginal ray in C */
  1442.  
  1443.        INTRIG_trace_line(3, clear_aperture / 2.0);
  1444.        od_cline = object_distance;
  1445.  
  1446.        /* Trace marginal ray in F */
  1447.  
  1448.        INTRIG_trace_line(6, clear_aperture / 2.0);
  1449.        od_fline = object_distance;
  1450.  
  1451.        aberr_lspher = od_sa[1][0] - od_sa[0][0];
  1452.        aberr_osc = 1.0 - (od_sa[1][0] * od_sa[1][1]) /
  1453.           (INTRIG_sin(od_sa[0][1]) * od_sa[0][0]);
  1454.        aberr_lchrom = od_fline - od_cline;
  1455.        max_lspher = INTRIG_sin(od_sa[0][1]);
  1456.  
  1457.        /* D light */
  1458.  
  1459.        max_lspher = 0.0000926 / (max_lspher * max_lspher);
  1460.        max_osc = 0.0025;
  1461.        max_lchrom = max_lspher;
  1462.     }
  1463.  
  1464. }
  1465.